MAPEO de COVID-19¶

Base de datos¶

In [1]:
import pandas as pd
df = pd.read_csv("C:/Users/Michael Encalada/Documents/GitHub/Magallanes/positivos_covid (1).csv", delimiter=';')
In [2]:
print(df.head())

   FECHA_CORTE DEPARTAMENTO  PROVINCIA     DISTRITO METODODX  EDAD       SEXO  \
0     20241203       TUMBES     TUMBES       TUMBES       AG  46.0   FEMENINO   
1     20241203         LIMA       LIMA  JESUS MARIA       AG  69.0   FEMENINO   
2     20241203   SAN MARTIN  MOYOBAMBA    MOYOBAMBA       AG  55.0   FEMENINO   
3     20241203     AREQUIPA   CAYLLOMA    COPORAQUE       AG  50.0  MASCULINO   
4     20241203         LIMA       LIMA  JESUS MARIA       AG  58.0  MASCULINO   

   FECHA_RESULTADO    UBIGEO  id_persona  
0       20221207.0  240101.0    203499.0  
1       20230822.0  150113.0    221397.0  
2       20240108.0  220101.0    295651.0  
3       20230824.0   40506.0    851625.0  
4       20221217.0  150113.0    287786.0
In [3]:
df=df[df.EDAD>50]
In [4]:
df.METODODX.value_counts()
Out[4]:
AG     599842
PCR    400134
PR     293382
Name: METODODX, dtype: int64
In [5]:
# NUEVA VARIABLE "AÑO"

df['FECHA_RESULTADO'] = pd.to_datetime(df['FECHA_RESULTADO'], format='%Y%m%d', errors='coerce')
In [6]:
df['AÑO'] = df['FECHA_RESULTADO'].dt.year
In [7]:
df.AÑO.head()
Out[7]:
1    2023.0
2    2024.0
4    2022.0
7    2023.0
9    2020.0
Name: AÑO, dtype: float64
In [8]:
df['AÑO'] = df['AÑO'].astype('Int64') #a numeros enteros
In [9]:
print(df[['AÑO']].head())

    AÑO
1  2023
2  2024
4  2022
7  2023
9  2020

Agreggate¶

In [10]:
# AG: Antigenos
# PCR: Prueba molecular
# PR: Prueba rápida serológica
In [11]:
df = df[df['AÑO'] != 1899]
In [12]:
df = df[df['PROVINCIA'] != 'EN INVESTIGACIÓN']
In [13]:
indexList=['AÑO','DEPARTAMENTO','PROVINCIA','METODODX']
aggregator={'METODODX':[len]}
prov=df.groupby(indexList,observed=True).agg(aggregator)
prov
Out[13]:
METODODX
len
AÑO DEPARTAMENTO PROVINCIA METODODX
2020 AMAZONAS BAGUA PCR 203
PR 2368
BONGARA PCR 20
PR 92
CHACHAPOYAS PCR 90
... ... ... ... ...
2024 TUMBES TUMBES PCR 15
ZARUMILLA PCR 6
UCAYALI CORONEL PORTILLO AG 16
PCR 17
PADRE ABAD AG 1

1886 rows × 1 columns

In [14]:
#wide
Draft=prov.unstack(3).fillna(0) #leftmost index in rows
Draft
Out[14]:
METODODX
len
METODODX AG PCR PR
AÑO DEPARTAMENTO PROVINCIA
2020 AMAZONAS BAGUA 0.0 203.0 2368.0
BONGARA 0.0 20.0 92.0
CHACHAPOYAS 0.0 90.0 419.0
CONDORCANQUI 0.0 14.0 540.0
LUYA 0.0 1.0 89.0
... ... ... ... ... ...
2024 TUMBES CONTRALMIRANTE VILLAR 0.0 4.0 0.0
TUMBES 11.0 15.0 0.0
ZARUMILLA 0.0 6.0 0.0
UCAYALI CORONEL PORTILLO 16.0 17.0 0.0
PADRE ABAD 1.0 0.0 0.0

920 rows × 3 columns

In [15]:
# % de personas con prueba AG
Draft['AG_pct']=Draft.iloc[:,1]/(Draft.iloc[:,0] + Draft.iloc[:,1])
prov_Ag=Draft['AG_pct'].unstack('AÑO').fillna(0)
prov_Ag
Out[15]:
AÑO 2020 2021 2022 2023 2024
DEPARTAMENTO PROVINCIA
AMAZONAS BAGUA 1.000000 0.144654 0.233813 0.240000 0.470588
BONGARA 1.000000 0.291971 0.466019 0.454545 0.928571
CHACHAPOYAS 1.000000 0.448919 0.414784 0.219178 0.322581
CONDORCANQUI 1.000000 0.033333 0.025000 0.000000 0.000000
LUYA 1.000000 0.036585 0.059361 0.111111 0.750000
... ... ... ... ... ... ...
TUMBES ZARUMILLA 0.923077 0.315556 0.492958 0.588235 1.000000
UCAYALI ATALAYA 0.000000 0.000000 0.000000 0.000000 0.000000
CORONEL PORTILLO 0.996875 0.186831 0.503120 0.545455 0.515152
PADRE ABAD 0.000000 0.078512 0.193694 0.000000 0.000000
PURUS 0.000000 0.000000 0.000000 1.000000 0.000000

196 rows × 5 columns

In [16]:
#data type
prov_Ag.columns #esta como numero los años
Out[16]:
Index([2020, 2021, 2022, 2023, 2024], dtype='Int64', name='AÑO')
In [17]:
prov_Ag.columns=['year'+str(x) for x in prov_Ag.columns]
In [18]:
prov_Ag
Out[18]:
year2020 year2021 year2022 year2023 year2024
DEPARTAMENTO PROVINCIA
AMAZONAS BAGUA 1.000000 0.144654 0.233813 0.240000 0.470588
BONGARA 1.000000 0.291971 0.466019 0.454545 0.928571
CHACHAPOYAS 1.000000 0.448919 0.414784 0.219178 0.322581
CONDORCANQUI 1.000000 0.033333 0.025000 0.000000 0.000000
LUYA 1.000000 0.036585 0.059361 0.111111 0.750000
... ... ... ... ... ... ...
TUMBES ZARUMILLA 0.923077 0.315556 0.492958 0.588235 1.000000
UCAYALI ATALAYA 0.000000 0.000000 0.000000 0.000000 0.000000
CORONEL PORTILLO 0.996875 0.186831 0.503120 0.545455 0.515152
PADRE ABAD 0.000000 0.078512 0.193694 0.000000 0.000000
PURUS 0.000000 0.000000 0.000000 1.000000 0.000000

196 rows × 5 columns

In [19]:
prov_Ag.reset_index(inplace=True)
prov_Ag
Out[19]:
DEPARTAMENTO PROVINCIA year2020 year2021 year2022 year2023 year2024
0 AMAZONAS BAGUA 1.000000 0.144654 0.233813 0.240000 0.470588
1 AMAZONAS BONGARA 1.000000 0.291971 0.466019 0.454545 0.928571
2 AMAZONAS CHACHAPOYAS 1.000000 0.448919 0.414784 0.219178 0.322581
3 AMAZONAS CONDORCANQUI 1.000000 0.033333 0.025000 0.000000 0.000000
4 AMAZONAS LUYA 1.000000 0.036585 0.059361 0.111111 0.750000
... ... ... ... ... ... ... ...
191 TUMBES ZARUMILLA 0.923077 0.315556 0.492958 0.588235 1.000000
192 UCAYALI ATALAYA 0.000000 0.000000 0.000000 0.000000 0.000000
193 UCAYALI CORONEL PORTILLO 0.996875 0.186831 0.503120 0.545455 0.515152
194 UCAYALI PADRE ABAD 0.000000 0.078512 0.193694 0.000000 0.000000
195 UCAYALI PURUS 0.000000 0.000000 0.000000 1.000000 0.000000

196 rows × 7 columns

MAPA¶

In [20]:
mapLink='https://github.com/SocialAnalytics-StrategicIntelligence/GeoDF_Analytics/raw/main/maps/ProvsINEI2023.zip'
#panadas vienen como zip hy, r no puede abrir por si acaso

import geopandas as gpd

provmap=gpd.read_file(mapLink)

provmap.info()

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 196 entries, 0 to 195
Data columns (total 6 columns):
 #   Column      Non-Null Count  Dtype   
---  ------      --------------  -----   
 0   OBJECTID    196 non-null    float64 
 1   CCDD        196 non-null    object  
 2   CCPP        196 non-null    object  
 3   DEPARTAMEN  196 non-null    object  
 4   PROVINCIA   196 non-null    object  
 5   geometry    196 non-null    geometry
dtypes: float64(1), geometry(1), object(4)
memory usage: 9.3+ KB
In [21]:
provmap['location']=['+'.join(x[0]) for x in zip(provmap.iloc[:,3:5].values)]
provmap.head(10)
Out[21]:
OBJECTID CCDD CCPP DEPARTAMEN PROVINCIA geometry location
0 1.0 01 01 AMAZONAS CHACHAPOYAS POLYGON ((-77.72614 -5.94354, -77.72486 -5.943... AMAZONAS+CHACHAPOYAS
1 2.0 01 02 AMAZONAS BAGUA POLYGON ((-78.61909 -4.51001, -78.61802 -4.510... AMAZONAS+BAGUA
2 3.0 01 03 AMAZONAS BONGARA POLYGON ((-77.72759 -5.14030, -77.72361 -5.140... AMAZONAS+BONGARA
3 4.0 01 04 AMAZONAS CONDORCANQUI POLYGON ((-77.81399 -2.99278, -77.81483 -2.995... AMAZONAS+CONDORCANQUI
4 5.0 01 05 AMAZONAS LUYA POLYGON ((-78.13023 -5.90370, -78.13011 -5.904... AMAZONAS+LUYA
5 6.0 01 06 AMAZONAS RODRIGUEZ DE MENDOZA POLYGON ((-77.44452 -6.05002, -77.44387 -6.050... AMAZONAS+RODRIGUEZ DE MENDOZA
6 7.0 01 07 AMAZONAS UTCUBAMBA POLYGON ((-78.09288 -5.36258, -78.09288 -5.364... AMAZONAS+UTCUBAMBA
7 8.0 02 01 ANCASH HUARAZ POLYGON ((-77.39870 -9.35563, -77.39852 -9.356... ANCASH+HUARAZ
8 9.0 02 02 ANCASH AIJA POLYGON ((-77.61368 -9.64900, -77.61241 -9.649... ANCASH+AIJA
9 10.0 02 03 ANCASH ANTONIO RAYMONDI POLYGON ((-77.08856 -8.97496, -77.08804 -8.975... ANCASH+ANTONIO RAYMONDI
In [22]:
print(prov_Ag.columns)

Index(['DEPARTAMENTO', 'PROVINCIA', 'year2020', 'year2021', 'year2022',
       'year2023', 'year2024'],
      dtype='object')
In [23]:
# Unir las columnas 'departamento' y 'provincia' en una nueva columna 'location'
prov_Ag['location'] = prov_Ag.apply(lambda row: f"{row['DEPARTAMENTO']} + {row['PROVINCIA']}", axis=1)

# Ver las primeras filas para confirmar
prov_Ag.head()
Out[23]:
DEPARTAMENTO PROVINCIA year2020 year2021 year2022 year2023 year2024 location
0 AMAZONAS BAGUA 1.0 0.144654 0.233813 0.240000 0.470588 AMAZONAS + BAGUA
1 AMAZONAS BONGARA 1.0 0.291971 0.466019 0.454545 0.928571 AMAZONAS + BONGARA
2 AMAZONAS CHACHAPOYAS 1.0 0.448919 0.414784 0.219178 0.322581 AMAZONAS + CHACHAPOYAS
3 AMAZONAS CONDORCANQUI 1.0 0.033333 0.025000 0.000000 0.000000 AMAZONAS + CONDORCANQUI
4 AMAZONAS LUYA 1.0 0.036585 0.059361 0.111111 0.750000 AMAZONAS + LUYA
In [24]:
## Limpieza
In [25]:
#!pip install unidecode
In [26]:
import unidecode

byePunctuation=lambda x: unidecode.unidecode(x)
prov_Ag['location']=prov_Ag['location'].apply(byePunctuation)
provmap['location']=provmap['location'].apply(byePunctuation)
In [27]:
prov_Ag['location']=prov_Ag.location.str.replace("\-|\_|\s+","",regex=True)
provmap['location']=provmap.location.str.replace("\-|\_|\s+","",regex=True)
In [28]:
nomatch_df=set(prov_Ag.location)- set(provmap.location)
nomatch_gdf=set(provmap.location)-set(prov_Ag.location)
In [29]:
len(nomatch_df), len(nomatch_gdf) #!!!
Out[29]:
(2, 2)
In [30]:
#!pip install thefuzz
In [31]:
# pick the closest match from nomatch_gdf for a value in nomatch_df
from thefuzz import process
[(dis,process.extractOne(dis,nomatch_gdf)) for dis in sorted(nomatch_df)]
Out[31]:
[('ANCASH+ANTONIORAIMONDI', ('ANCASH+ANTONIORAYMONDI', 95)),
 ('ICA+NAZCA', ('ICA+NASCA', 89))]
In [32]:
# Diccionario de cambios
{dis:process.extractOne(dis,nomatch_gdf)[0] for dis in sorted(nomatch_df)}
Out[32]:
{'ANCASH+ANTONIORAIMONDI': 'ANCASH+ANTONIORAYMONDI', 'ICA+NAZCA': 'ICA+NASCA'}
In [33]:
# then:
changesinDF={dis:process.extractOne(dis,nomatch_gdf)[0] for dis in sorted(nomatch_df)}
In [34]:
prov_Ag.replace({'location': changesinDF}, inplace=True)
In [35]:
# comprobamos
In [36]:
nomatch_df=set(prov_Ag.location)- set(provmap.location)
nomatch_gdf=set(provmap.location)-set(prov_Ag.location)

[(dis,process.extractOne(dis,nomatch_gdf)) for dis in sorted(nomatch_df)]
Out[36]:
[]
In [37]:
#merge
prov_Ag_map=provmap.merge(prov_Ag, on='location',how='left',indicator='flag')
In [38]:
prov_Ag_map.info()

<class 'geopandas.geodataframe.GeoDataFrame'>
Int64Index: 196 entries, 0 to 195
Data columns (total 15 columns):
 #   Column        Non-Null Count  Dtype   
---  ------        --------------  -----   
 0   OBJECTID      196 non-null    float64 
 1   CCDD          196 non-null    object  
 2   CCPP          196 non-null    object  
 3   DEPARTAMEN    196 non-null    object  
 4   PROVINCIA_x   196 non-null    object  
 5   geometry      196 non-null    geometry
 6   location      196 non-null    object  
 7   DEPARTAMENTO  196 non-null    object  
 8   PROVINCIA_y   196 non-null    object  
 9   year2020      196 non-null    float64 
 10  year2021      196 non-null    float64 
 11  year2022      196 non-null    float64 
 12  year2023      196 non-null    float64 
 13  year2024      196 non-null    float64 
 14  flag          196 non-null    category
dtypes: category(1), float64(6), geometry(1), object(7)
memory usage: 23.3+ KB
In [39]:
prov_Ag_map['flag']=prov_Ag_map.flag.astype(str)
In [40]:
#
bye=['DEPARTAMEN', 'PROVINCIA_x', 'CCPP','CCDD']
prov_Ag_map.drop(columns=bye,inplace=True)
prov_Ag_map.head()
Out[40]:
OBJECTID geometry location DEPARTAMENTO PROVINCIA_y year2020 year2021 year2022 year2023 year2024 flag
0 1.0 POLYGON ((-77.72614 -5.94354, -77.72486 -5.943... AMAZONAS+CHACHAPOYAS AMAZONAS CHACHAPOYAS 1.0 0.448919 0.414784 0.219178 0.322581 both
1 2.0 POLYGON ((-78.61909 -4.51001, -78.61802 -4.510... AMAZONAS+BAGUA AMAZONAS BAGUA 1.0 0.144654 0.233813 0.240000 0.470588 both
2 3.0 POLYGON ((-77.72759 -5.14030, -77.72361 -5.140... AMAZONAS+BONGARA AMAZONAS BONGARA 1.0 0.291971 0.466019 0.454545 0.928571 both
3 4.0 POLYGON ((-77.81399 -2.99278, -77.81483 -2.995... AMAZONAS+CONDORCANQUI AMAZONAS CONDORCANQUI 1.0 0.033333 0.025000 0.000000 0.000000 both
4 5.0 POLYGON ((-78.13023 -5.90370, -78.13011 -5.904... AMAZONAS+LUYA AMAZONAS LUYA 1.0 0.036585 0.059361 0.111111 0.750000 both
In [41]:
# lidiar con los CEROS
prov_Ag_map.fillna(0,inplace=True)

Guardar el geoDF

In [42]:
#crear carpeta maps, donde guardar

import os
prov_Ag_map.to_file(os.path.join('C:/Users/Michael Encalada/Documents/GitHub/Magallanes/week10_spatial/maps/',"provinciasPeru.gpkg"), layer='provinciasDengue', driver="GPKG")

Explorar variables¶

In [43]:
# statistics
prov_Ag_map.year2022.describe()
Out[43]:
count    196.000000
mean       0.125021
std        0.161102
min        0.000000
25%        0.016327
50%        0.068471
75%        0.170919
max        1.000000
Name: year2022, dtype: float64
In [44]:
#! pip install seaborn
In [45]:
# grafico
import seaborn as sea

sea.boxplot(prov_Ag_map.year2022, color='yellow',orient='h')
Out[45]:
<Axes: xlabel='year2022'>
No description has been provided for this image

Imterpretación: El 74% de las provincias no llegan al 40% de realización de pruebas antigenas para detectar covid 19.

In [46]:
#grafica de trasnformacion cuantilica?
from sklearn.preprocessing import QuantileTransformer
qt = QuantileTransformer(n_quantiles=100, random_state=0,output_distribution='normal')
qt_result=qt.fit_transform(prov_Ag_map[['year2022']])
sea.boxplot(qt_result, color='yellow',orient='h')
Out[46]:
<Axes: >
No description has been provided for this image

Explicación: Transofrmador quantilico: evita que los atipicos sesguen los resultados No estas eliminando los atipicos, lo estas normalizando haciendolo parte de la distribución. En todo la politica publica, hay que pensar en lo intervinible. Necesito algo que sea operable. Lo hago viendo al grupo donde se puede trabajr. Eso es el amarillo

In [47]:
prov_Ag_map['year_2022_qt']= qt_result

Correlación Espacial¶

Vecindario

In [48]:
#!pip install pysal 
#conda install -c conda-forge pysal
In [49]:
from libpysal.weights import Queen, Rook, KNN
In [50]:
# rook
w_rook = Rook.from_dataframe(prov_Ag_map,use_index=False)
In [51]:
# queen
w_queen = Queen.from_dataframe(prov_Ag_map,use_index=False)
In [52]:
# k nearest neighbors
w_knn = KNN.from_dataframe(prov_Ag_map, k=8) #quiero 8 mas cercanos a mi.

Moran's correlation¶

In [53]:
# needed for spatial correlation
w_queen.transform = 'R'
In [54]:
pd.DataFrame(*w_queen.full()).sum(axis=1) # 1 means both are neighbors
Out[54]:
0      1.0
1      1.0
2      1.0
3      1.0
4      1.0
      ... 
191    1.0
192    1.0
193    1.0
194    1.0
195    1.0
Length: 196, dtype: float64

INDICE DE MORAN

Quiero ver si hay correlacion entre mi vecindario con dengue. Caunto dengue hay en mi barro

ojo: que a veces la frecuencia entre personas (redes) va mas alla del espacio. frecuentas con un grupo a pesar de distancias.

In [55]:
from esda.moran import Moran

moranCOVID = Moran(prov_Ag_map['year_2022_qt'], w_queen)
moranCOVID.I,moranCOVID.p_sim

#correlacion a 0.51 y es significativa por es 0.001 ***
Out[55]:
(0.19136370498258676, 0.001)

interpretacion: Si sale no significativo, significa que el virus no tiene una variable espacial. no tiene efecto La cercania no tiene una relacion con el territorio. Vota en contra de keiko, porque done vivo mis vecinos tambien votan en contra de keiko

Grafico

In [56]:
from splot.esda import moran_scatterplot
import matplotlib.pyplot as plt

fig, ax = moran_scatterplot(moranCOVID)
ax.set_xlabel('Covid_prueba_share')
ax.set_ylabel('SpatialLag_Covid_prueba_share')
Out[56]:
Text(0, 0.5, 'SpatialLag_Covid_prueba_share')
No description has been provided for this image

Interpretacion: a parte de la linea, si no tambien los cuadrantes Los que estan a la izquierda, abajo, son provincias donde el dengue esta bajo y hay cercania entre provincias bajo.

A la izquierda arriba, es dengue alto con vecinos cercanas con dengue alto. A la derecha abajo, no contagiados, casi asilados.

izquierda abajo COLD SPOT ARRIB derecha, hot spot

Otro grafico

In [59]:
# The scatterplot with local info

from esda.moran import Moran_Local

# calculate Moran_Local and plot
lisaCOVID = Moran_Local(y=prov_Ag_map['year_2022_qt'], w=w_knn,seed=2022)
fig, ax = moran_scatterplot(lisaCOVID,p=0.05)
ax.set_xlabel('Covid_prueba_share')
ax.set_ylabel('SpatialLag_Covid_prueba_share');
No description has been provided for this image

Mapa de peru

In [58]:
# the map with the spots and outliers

from splot.esda import lisa_cluster
f, ax = plt.subplots(1, figsize=(12, 12))
plt.title('Spots and Outliers')
fig = lisa_cluster(lisaCOVID,
                   prov_Ag_map,ax=ax,
                   legend_kwds={'loc': 'center left',
                                'bbox_to_anchor': (0.7, 0.6)})
No description has been provided for this image

Agregar data to my gdf

In [60]:
# quadrant
lisaCOVID.q
Out[60]:
array([1, 1, 1, 1, 1, 1, 1, 4, 3, 3, 3, 3, 3, 3, 4, 3, 3, 4, 4, 3, 3, 3,
       3, 3, 4, 3, 3, 1, 4, 1, 4, 1, 4, 1, 4, 3, 3, 3, 4, 3, 4, 3, 1, 3,
       3, 1, 1, 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 3, 2, 2, 1, 2, 2, 1, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 4, 3, 3, 4, 3, 1, 1, 4, 3, 1, 1, 4, 4, 1, 1, 1, 1, 1, 2, 1, 1,
       2, 1, 1, 1, 1, 2, 4, 1, 1, 3, 4, 4, 1, 4, 1, 1, 1, 1, 4, 3, 1, 1,
       1, 1, 4, 3, 1, 1, 1, 1, 4, 1, 2, 1, 2, 4, 2, 4, 4, 2, 4, 1, 4, 1,
       4, 1, 2, 4, 1, 1, 1, 2, 4, 3, 2, 4, 4, 4, 4, 1, 2, 3, 4, 4, 4, 1,
       4, 2, 1, 1, 1, 1, 1, 1, 4, 4, 3, 3, 3, 1, 1, 1, 4, 2, 1, 3])
In [61]:
# significance
lisaCOVID.p_sim
Out[61]:
array([0.165, 0.301, 0.099, 0.016, 0.051, 0.129, 0.029, 0.209, 0.087,
       0.002, 0.012, 0.012, 0.013, 0.004, 0.335, 0.11 , 0.005, 0.09 ,
       0.017, 0.013, 0.079, 0.13 , 0.002, 0.017, 0.132, 0.003, 0.02 ,
       0.152, 0.181, 0.441, 0.058, 0.114, 0.239, 0.075, 0.404, 0.084,
       0.019, 0.098, 0.181, 0.066, 0.157, 0.132, 0.461, 0.333, 0.179,
       0.023, 0.332, 0.009, 0.02 , 0.008, 0.031, 0.052, 0.164, 0.135,
       0.028, 0.052, 0.005, 0.121, 0.005, 0.101, 0.083, 0.142, 0.066,
       0.049, 0.047, 0.025, 0.04 , 0.298, 0.22 , 0.225, 0.451, 0.232,
       0.301, 0.269, 0.26 , 0.265, 0.243, 0.38 , 0.305, 0.195, 0.021,
       0.039, 0.029, 0.095, 0.014, 0.469, 0.015, 0.406, 0.5  , 0.059,
       0.001, 0.013, 0.077, 0.006, 0.498, 0.015, 0.073, 0.408, 0.499,
       0.016, 0.37 , 0.331, 0.015, 0.036, 0.173, 0.214, 0.172, 0.217,
       0.028, 0.271, 0.223, 0.095, 0.085, 0.184, 0.035, 0.014, 0.373,
       0.036, 0.145, 0.028, 0.231, 0.2  , 0.066, 0.309, 0.03 , 0.067,
       0.333, 0.036, 0.017, 0.091, 0.242, 0.103, 0.475, 0.016, 0.159,
       0.034, 0.017, 0.307, 0.144, 0.056, 0.369, 0.264, 0.093, 0.046,
       0.031, 0.358, 0.11 , 0.141, 0.039, 0.357, 0.058, 0.264, 0.056,
       0.112, 0.31 , 0.15 , 0.001, 0.234, 0.159, 0.003, 0.006, 0.122,
       0.03 , 0.332, 0.456, 0.012, 0.081, 0.41 , 0.369, 0.481, 0.275,
       0.327, 0.152, 0.392, 0.044, 0.083, 0.034, 0.136, 0.137, 0.21 ,
       0.465, 0.42 , 0.078, 0.438, 0.001, 0.036, 0.218, 0.069, 0.149,
       0.009, 0.011, 0.011, 0.327, 0.488, 0.333, 0.207])
In [63]:
# quadrant: 1 HH,  2 LH,  3 LL,  4 HL
pd.Series(lisaCOVID.q).value_counts()
Out[63]:
1    94
4    43
3    42
2    17
dtype: int64

Agregar a la correlación espacial

In [64]:
prov_Ag_map['COVID_quadrant']=[l if p <0.05 else 0 for l,p in zip(lisaCOVID.q,lisaCOVID.p_sim)  ]
prov_Ag_map['COVID_quadrant'].value_counts()
Out[64]:
0    131
1     32
3     21
4      9
2      3
Name: COVID_quadrant, dtype: int64
In [66]:
#recodificar
labels = [ '0 no_sig', '1 hotSpot', '2 coldOutlier', '3 coldSpot', '4 hotOutlier']

prov_Ag_map['COVID_quadrant_names']=[labels[i] for i in prov_Ag_map['COVID_quadrant']]

prov_Ag_map['COVID_quadrant_names'].value_counts()
Out[66]:
0 no_sig         131
1 hotSpot         32
3 coldSpot        21
4 hotOutlier       9
2 coldOutlier      3
Name: COVID_quadrant_names, dtype: int64

Regraficar

In [67]:
from matplotlib import colors
myColMap = colors.ListedColormap([ 'ghostwhite', 'red', 'green', 'black','orange'])

f, ax = plt.subplots(1, figsize=(12,12))


plt.title('Spots and Outliers')

prov_Ag_map.plot(column='COVID_quadrant_names', 
                categorical=True,
                cmap=myColMap,
                linewidth=0.1, 
                edgecolor='white',
                legend=True,
                legend_kwds={'loc': 'center left', 
                             'bbox_to_anchor': (0.7, 0.6)},
                ax=ax)
# Remove axis
ax.set_axis_off()
# Display the map
plt.show()
No description has been provided for this image
In [68]:
prov_Ag_map.explore("COVID_quadrant_names", categorical=True,tooltip='location',cmap=myColMap)
Out[68]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [70]:
import folium

map1=prov_Ag_map[prov_Ag_map.COVID_quadrant_names=='1 hotSpot']
map2=prov_Ag_map[prov_Ag_map.COVID_quadrant_names=='2 coldOutlier']
map3=prov_Ag_map[prov_Ag_map.COVID_quadrant_names=='3 coldSpot']
map4=prov_Ag_map[prov_Ag_map.COVID_quadrant_names=='4 hotOutlier']

m = map1.explore(
    color="red",  
    tooltip=False,  # hide tooltip
    popup=["location"],  # (on-click)
    name="hotSpot"  # name of the layer in the map
)

map2.explore(
    m=m, # notice
    color="green",  
    tooltip=False,  
    popup=["location"],
    name="coldOutlier"
)

map3.explore(
    m=m,
    color="black",  
    tooltip=False,  
    popup=["location"],
    name="coldSpot", 
)

map4.explore(
    m=m,
    color="orange", 
    tooltip=False,  
    popup=["location"],
    name="hotOutlier", 
)

folium.TileLayer("CartoDB positron", show=False).add_to(m)  # use folium to add alternative tiles
folium.LayerControl(collapsed=True).add_to(m)  # use folium to add layer control

m  # show map
Out[70]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Bivariate LISA¶

In [71]:
#from esda.moran import Moran_BV, Moran_Local_BV
from esda.moran import Moran_BV

mbi = Moran_BV(prov_Ag_map['year2021'],  prov_Ag_map['year2022'],  w_queen)
mbi.I,mbi.p_sim
Out[71]:
(0.14224711545767468, 0.003)
In [72]:
# The scatterplot with local info
from esda.moran import Moran_Local_BV

# calculate Moran_Local and plot
lisaCOVID_bv = Moran_Local_BV(y=prov_Ag_map['year2021'],
                               x=prov_Ag_map['year2022'],
                               w=w_queen)

fig, ax = moran_scatterplot(lisaCOVID_bv, p=0.05,aspect_equal=True)

ax.set_xlabel('COVID_2022')
ax.set_ylabel('SpatialLag_COVID_2021')
plt.show()
No description has been provided for this image
In [73]:
prov_Ag_map['COVID_quadrant_21_22']=[l if p <0.05 else 0 for l,p in zip(lisaCOVID_bv.q,lisaCOVID_bv.p_sim)  ]

labels = [ '0 no_sig', '1 hotSpot', '2 coldOutlier', '3 coldSpot', '4 hotOutlier']

prov_Ag_map['COVID_quadrant_21_22_names']=[labels[i] for i in prov_Ag_map['COVID_quadrant_21_22']]
In [74]:
# see new columns
prov_Ag_map
Out[74]:
OBJECTID geometry location DEPARTAMENTO PROVINCIA_y year2020 year2021 year2022 year2023 year2024 flag year_2022_qt COVID_quadrant COVID_quadrant_names COVID_quadrant_21_22 COVID_quadrant_21_22_names
0 1.0 POLYGON ((-77.72614 -5.94354, -77.72486 -5.943... AMAZONAS+CHACHAPOYAS AMAZONAS CHACHAPOYAS 1.000000 0.448919 0.414784 0.219178 0.322581 both 1.478159 0 0 no_sig 0 0 no_sig
1 2.0 POLYGON ((-78.61909 -4.51001, -78.61802 -4.510... AMAZONAS+BAGUA AMAZONAS BAGUA 1.000000 0.144654 0.233813 0.240000 0.470588 both 1.027998 0 0 no_sig 0 0 no_sig
2 3.0 POLYGON ((-77.72759 -5.14030, -77.72361 -5.140... AMAZONAS+BONGARA AMAZONAS BONGARA 1.000000 0.291971 0.466019 0.454545 0.928571 both 1.720078 0 0 no_sig 0 0 no_sig
3 4.0 POLYGON ((-77.81399 -2.99278, -77.81483 -2.995... AMAZONAS+CONDORCANQUI AMAZONAS CONDORCANQUI 1.000000 0.033333 0.025000 0.000000 0.000000 both -0.502255 1 1 hotSpot 0 0 no_sig
4 5.0 POLYGON ((-78.13023 -5.90370, -78.13011 -5.904... AMAZONAS+LUYA AMAZONAS LUYA 1.000000 0.036585 0.059361 0.111111 0.750000 both -0.140240 0 0 no_sig 2 2 coldOutlier
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
191 192.0 POLYGON ((-80.28521 -3.41276, -80.28406 -3.412... TUMBES+ZARUMILLA TUMBES ZARUMILLA 0.923077 0.315556 0.492958 0.588235 1.000000 both 1.827310 1 1 hotSpot 0 0 no_sig
192 193.0 POLYGON ((-74.47145 -7.27617, -74.47052 -7.277... UCAYALI+CORONELPORTILLO UCAYALI CORONEL PORTILLO 0.996875 0.186831 0.503120 0.545455 0.515152 both 1.866251 0 0 no_sig 4 4 hotOutlier
193 194.0 POLYGON ((-73.18146 -9.41174, -73.13475 -9.411... UCAYALI+ATALAYA UCAYALI ATALAYA 0.000000 0.000000 0.000000 0.000000 0.000000 both -5.199338 0 0 no_sig 0 0 no_sig
194 195.0 POLYGON ((-75.43663 -8.22999, -75.43651 -8.230... UCAYALI+PADREABAD UCAYALI PADRE ABAD 0.000000 0.078512 0.193694 0.000000 0.000000 both 0.778268 0 0 no_sig 0 0 no_sig
195 196.0 POLYGON ((-70.61380 -9.87339, -70.62140 -9.878... UCAYALI+PURUS UCAYALI PURUS 0.000000 0.000000 0.000000 1.000000 0.000000 both -5.199338 0 0 no_sig 0 0 no_sig

196 rows × 16 columns

Un ultimo mapita

In [75]:
from matplotlib import colors
myColMap = colors.ListedColormap([ 'ghostwhite', 'red', 'green', 'black','orange'])


f, ax = plt.subplots(1, figsize=(12,12))

plt.title('Spots and Outliers')

prov_Ag_map.plot(column='COVID_quadrant_21_22_names', 
                categorical=True,
                cmap=myColMap,
                linewidth=0.1, 
                edgecolor='white',
                legend=True,
                legend_kwds={'loc': 'center left', 
                             'bbox_to_anchor': (0.7, 0.6)},
                ax=ax)
# Remove axis
ax.set_axis_off()
# Display the map
plt.show()
No description has been provided for this image

Convertir a Html¶

FIN¶